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ABSTRACT 

A new Gauss-Newton algorithm is presented for solving non- 
linear least squares problems. The problem statement may 
include simple bounds or more general constraints on the 
unknowns. The algorithm uses a trust region that allows the 
objective function to increase with logic for retreating to 
best values. The computations for the linear problem are 
done using a least squares system solver that allows for 
simple bounds and linear constraints. The trust region 
limits are defined by a box around the current point. In 
its current form the algorithm is effective only for prob- 
lems with small residuals# linear constraints and dense 
Jacobian matrices. Results on a set of test problems are 
encouraging. 

t The research of this author was carried out at the Jet Pro- 
pulsion Laboratory, California Institute of Technology# 
under contract with the National Aeronautics and Space 
" dministration . 
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Introduction 

The constrained nonlinear least squares problem which we consider is 
stated as follows. 

Problem CNLLS 

Let f (x) be a differentiable mapping of a rectangular 
domain T = {x:aj £ Xj £ Bjr j = lf.../N} into R^, and g(x) 
(1) be a differentiable mapping of T into R^'. Find an x e T, 

' A ' ”2 

“i £ 9i^5) £ ®i» i-l».*.rl*» such that llf(x)ll2 = ^ fi(x) 
is minimized at x = x. 


In the statement of problem CNLLS we emphasize that the functions 
f(x) and g (x) are defined and differentiable for every xeT. 

Problems for which f(x) is defined only when xeT and a' < g(x) < fi* 

are not included. The solution x will satisfy the bounds on g(x), if 

possibler but g(x) may be evaluated at an arbitrary point xeT. 


The following notation is used, 


?p 

X on the p-th iteration 


f(x) on the p-th iteration 

9p 

g(x) on the p-th iteration 

Gp 

3g/3x on the p-th iteration 

Jp 

3f/3x on the p-th iteration 

Axp 

Xp - Xp+i 


Jhe basic iteration, known as the (constrained) Gauss-Newton method, 
is given by the sequence of constrained linear least squat es problems 





Probably the most challenging aspect of developing an algorithm, 
based on (3), is in making the method converge with good efficiency. 
(Our notion of efficiency is based on the usual one, the number of 
function and derivative evaluations, fp, Jp, gp, and Gp . The cost 
of the numerical linear algebra to solve the constrained linear least 
squares problem of (3) is assumed to be less important than the cost 
of evaluating these functions and derivr'tives. ) It is frequently 
true that the iteration (3) will fail to converge if HAxpli 
is not restricted in some way. Too tight a restriction, however, 
will decrease the efficiency. 

We describe an algorithm for bounding liAxpll which (based on 
performance) gives better efficiency, on the average, than the 
MINPACK algorithm, [1]. Cur work has not yet resulted in a finished 
software product. What we are reporting here is the fact that our 
basic methods look promising. 

It must be emphasized that even when Xp is close to a solution, 

the iteration (3) may be unacceptably slow or may not converge. 
Cases A-D below must be separately considered. They can.- of course, 
appear in any combination. We feel that our algorithm currently 

performs reasonably in Case A. 

A. The matrix Jp is full rank and iifpi!-^c-0. Both Jjfp and Axp con- 
verge quadra^ically to zero. (A simple model is f(x)=sin x, 
near x=0. ) 

B. The matrix Jp is of full rank and lifpli-^ci >0. Iteration (3) gives 
only linear convergence to the solution. (A simple model 
is fT(x)=(l+sin X, 2-e^) , near x=0.) 

C. At the solution Jp is not of full rank, but this does not 

limit the amount nf(x)ii can be reduced. Iteration (3) 

converges linearly to the solution. (A simple model is f(x) = 
x sin x, near x=0.) 


D. At the solution Jp is not of full rank and this limits the 
amount itf(x)H can be reduced. Iteration (3) does not con- 

- A# M 

verge. (A simple model is f(x) « 1 t sin^x, near x«0.) 

We anticipate that in Part II of this report (not yet written) modi- 
fications of the iteration (3) will be presented. These modifica- 
tions will involve using (3) where it is effective and higher order 
terms where it is not. We believe that oi.r modified algorithm will 
be efficient for all the cases except for uase C when case D is also 
present. We believe this combination is a rare event in the computa- 
tions that will apply our methods. 

Our current algorithm is effective only on problems where Gp is a 
constant matrix. We expect to modify the algorithm for nonlinear 
constraints in Part II of this report. 

Section 2 gives a broad outline of the algorithm and contrasts it to 
the Levenberg-Marquardt method used, for example in MINPACK [1] and 
the package NL2SOL [2]. Section 3 gives the details of our algorithm 
and the rationale for various aspects of its design. Section 4 sum- 
marizes our results on the test problems and contains some discussion 
of the test set itself. 

2. Outline of the Algorithm 

in this section we summarize our approach for solving the non- 
linear least squares problem of (1) . The key to efficiency 
is in accepting steps Axp for which !tfp+ih > iifpii. Large moves 
that increase iif(x)ii are frequently good moves in the sense 

that they get closer to the solution. in turn, this requires 
logic so that the algorithm can retreat to its previous best 
values if it senses that making these moves got it into trou- 
ble. We provide logic for sensing trouble and for making the 
retreat. 

The algorithm begins with a user-given starting point, x^. An 

N-dimensional box, B, is constructed that contains x_. The 

■ -.p 

lengths of each dimension of B are chosen so that B remains in T 
and liAxpii remains "reasonable.” We will elaborate on 
this below and in the next section. 

There are five cases for altering these dimensions that we will 
discuss in detail in Section 3. These cases are 

1. Selecting the dimensions of B initially. 

2. Selecting the dimensions after successive best values of 
iif(x)ii are obtained. 

3. Selecting the dimensions after finding a best value of 
iif(x)ii after a value that is not a best. 


4. Selecting the dimensions after giving up and retreating to a 
best value for Hf(x)lt. 

5. Selecting the dimensions after obtaining a value for 
ljf(x}ii larger than the best obtained previously. 

The Levenberg-Marquardt method for soiiring nonlinear least squares 
problems (with no constraints) can be regarded as placing an ellip- 
soidal trust region about each point Xp. The MINPACK implementa- 
tion [1] of this method chooses the relative sizes of the axes of 
this ellipsoid based on the column norms for Jp. This type of 
scaling is appropriate when the error in each variable contributes 
equally co the residual » a condition that is frequently not 
satisified. 

Our algorithm uses adaptive scaling that depends on changes in the 
individual components of x. Namely, the relative weights used to 

compute the future moves is approximated by the size of moves made 
in the past. We tried column scaling using reciprocals of column 
norms at successive best x values. Based on overall performance, we 

prefer the adaptive scaling outlined above. A Levenberg-Marquardt 
algorithm that uses our approach for column scaling could be 
developed. 

Briefly stated, our decision to retreat to a best value for 
ilf(x)li is made in one of two ways: 

1. If ''JpiW4'"!<*BEST>" * “bESt'^*BEST-!bESt"> ' 

retreat. If this condition is not satisfied, the linear approx- 
imation is telling us that we still have reasonable expectation 
for improvement, despite the fact that Ilf (x) II is large. 

A, 

2. If llf(x)ll gets too large after a move from a best x value, 

the initial move was probably too large. We retreat if 

lifpil>f, where "f=2llfpll at the first step with llfpll^llfp_xH. This 

allows an arbitrarily large increase on the first step; after 
that llf(x)ll is bounded. The final version of our algorithm 

A, ^ 

is not expected to require this condition. 

The characteristics that we consider vital here are: An arbitrar ily 

large increase in Ilf (x) II is tolerated; there is no ad hoc limit 

on how many iterations will be taken before retreating; and retreat- 
ing immediately if the residual to the current linear problem does 
not hold the prospect of an eventual new best. 


The nonlinear l^att Bquarta tlqortthm using simple boun<tg to 
define the trust region . ' , ^ 

Oar trust region aaounte to a ehoiee forf the dinensions of the 

box that we place around x^. The dimensions of the box are 

~P 

varied, roughly speaking, so that when the function decreases 
the dimensions increase. Otherwise the dimensions tend to 
decrease. 

We introduce the following list of variables that are used to 
describe the details of the algorithm. 

X t.^e vector of independent variables at Step p, i.e., Xn. 

Xj^ component i of x, i«l,...,N 

6 i change on the current iteration, i.e., component i of 

Axp. 

f function whose Euclidean length (least squares norm) is 

being minimized. 

xq value of x from among all the previous Xp that gave the 
minimum of Ilf (x) II; i.e., the "best” x. 

fg smallest value of !lf(x)ll so far, f^o^lif (X3) II 

% value of llf(Xp_2)ll, last value for llf(x)li 

f(2 value of Hf(Xp)ii, current value for iif(x)ii 

b£ bound used in computing 62 . 

bound used in computing 5 ^ after a new best 

k count of iterations since a new best. (Set to 0 to flag 

initialization.) 

KASE if KASE >0 got a new best after KASE iterations. (Set to 

0 for initialization; set to -1 when we retreat to best x.) 

Ci length of column i of the Jacobian matrix Jp. 

a a scalar used in computing b^ 

P predicted value of fc based on solving the constrained 

linear problem 

boost, a boost factor for Increasing the bounds on suc- 
cessive best x's. 


algorithm. The language we 
use here is close to SFTRAKSf {3]. We first describe the algorithm 
and follow this by some details about why we chosa the specific fac- 
tors to vary the trust region. The reader is warned that some 
details of the algor tthiB given here are obsolete even at this 
writing. The results shown in Tables 1-4 do correspond to the 
following algorithm, however. 

H-K-1 Algorithm for Nonlinear Least Squares 

6 £ ® 0 , i*l , . . • ,N 

k = 0 

fg » 00 

P * 0 

retreat « .false. 

terminate = .false. 

DO FOREVER 

IF (retreat) THEN 

C We must retreat back to best x. 


k =0 
KASE = -1 

X » xb 

Else 

KASE = k 

Xj_ — X^ — 1*1 , . . • ,N 

If (terminate) Exit from Algorithm 
End if 

Compute J, f, G, and g at x. 

fc = ilf(x)il 

Test for convergence. 

If (terminate) Exit from Algorithm 

If (fc < fs) k = 0 

If (k » 0) then 

C We want to position at best x values. 


fB * fc 


DO CASE (2«KASB» 3} 


Case 1 

C We liniiied lately got a new best x. 

C The following formula for a has the following effect. 

C If 1 P' ct - 1 + 1/AI.FAC. 

C If f2 « P * fL» then a s 1. 

C If fg>>P * fL» then a » 0.5. 

a = max(fg - p2,0) 

a + P * (ALFAC + 1) * (f^ - P) 

“ “ 2a + P * ALFAC * (f £ y) 

ALFAC = 0.125 * (- + ALFAC) 

a 

boost * min ( (a+1) ^boostf 10^®) 

Case 2 

C We are at the initial x. 

ALFAC * 1/128 
Do For i « 1, . . . ,N 

Bi » -Xi 

If (Bi » 0) then 

If (Ci / 0) then 
Bi » “fc/^i 

Else 

Bi * -1 
End If 
End If 
bi « Bi 
End For 
XB * X 

a ■ 1 

boost » 0.5 

Exit from the "If (k « 0)" block. 

Case 3 

C We are retreating to best x. 


a •0.125 
boost ■0.25 

Case 

Hot linmed lately a best Je. 

a ■ 0.25 
boost ■ 1 

End Case 

Do For i “ If. 

«1 = f*B)i - Xi 

If ( 6 1 « 0) then 
bi ■ 

Else 

bj « a*SIGNC6i-tSi,5i) + boost* 6 

{xb) 1 » *i 
End If 

End For 
Else 

The nonlinear residual nc-ni is not a new best 
ALFAC -0.25 

If(k » 1) then 
f * 2*fc 

P « (fB+P)/2 
Else If (fQ>f) then 

The function norm is too big. 

retreat « .true. 

cycle on "Do Forever" loop 

End if 

End if 

a - (0.25*fc + 0.5*fi,)/(fc + ^L> 

Do For i « 1, . . . ,N 

If » 0) then 
bi » a*bj 



FJaT'^'FriHW! 


Compute the boXf with each dimension deilnCKl hy 

(ii»ui) that will be trust region ior t¥m 

problem. 


If (bl<d) then 
li - bi 

»1 - -£i* 


c 

c 


c 

c 


Restrict the trust region tounds further if the user has 
given bounds on the variables. 

x’i • max {li, xi - ei) 
ui » min (ui, xi - at) 

End For 

If the user has given general constraints , compute the 
bounds (£j, Uj) for the constraint equations. 

Do For j “ 1 r . . . , L 


C 

C 

C 


9j(x) " 

Jj ■ 9j(x) - Qtj 
End For 

Solve linear least squares prcolem 14] with bounds and 
generals constraints to get j ■ l....,N, and the 

linear residual, vector length, P. 



o u 


If (P > P) then 
retreat « .true. 

End if 

If (.not. retreat) th^n 

If (liAx il is small and certain other 
P 

conditions are satisfied) then 
terminate = .true. 

End if 
End if 

= *C 

Reduce bounds if the size of the move was well below what 
was allowed 


a = 0 

Do For j = 1, . . . ,N 

t = 6Vb, 

3 j 

If (t > 0) then 


a = max (a,t) 


Else 


a = max (a,-(9/8)*t) 


End If 


End For 


If (a £ 0.5) then 


0 = 2*0 


Do For j = 1» . . . ,N 


bj = 0 *bj 


If (k = 1) Bj = bj 


End For 


boost = 1 


End If 


End Forever 


14 


Jam'll.. "Njk'ii!' 'iiii iiL Jl lik 


He close this section with sone discussion of the reasons for our 
choice of the dinensions of the box» B. These are cases 1-5 given 
in Section 2. The algorithm tends to select slightly smaller bounds 
for those components which are changing direction. The direction is 
set at the initial point as the sign of the value. For example, a 
variable with the value 1 can move to points in the interval 
[1/9, 21. 

1. selecting the J3lmensiona of B Initially 

If a variable is nonzero, we assume that its value is in the 

approximate scale of a solution for this parameter. The user 

should avoid nonzero initial guesses if this is not the case. 
The dimension of the box, in this coordinate, is 17/9 the magni- 
tude of the variable. If the variable is zero and c£ 0, 
we use an approximation to a single variable Gauss-Newton step. 
This approximation assumes that the residual could be reduced to 
zero with that component. This yields (17/9) iifQii/c^ as 

the dimension of B. If » 0, we use the value of (17/9) for 

the dimension in this coordinate. These bounds tend to be large, 
especially if one were to insist on a reduction in tif(x}ii. 

Ne have not done much work on this choice of initial bounds. 
Perhaps the choice should be symmetric in a relative sense for 
example. More testing of this could be done, and it is likely 
that tuning would yield some improvements in the algorithm's 
performance. The choice we made was guided by retaining simplic- 
ity of implementation. If a choice of a nonzero initial value 
is off scale y overal orders of magnitude, our algorithm may 
be inefficient. We plan to make changes to the algorithm to 
deal with this out-of-scale problem at a later time. 

2. Selectin g th e Dimensions After Successiy e Best v^^ues gf 
lif(x)ii are obtained 


Our principle for choosing the scale on the size of B is based 
on past behavior. Specifically, the bounds are closely related 
to the distance between adjacent values of x. There are three 

further effects that we use to choose the dimensions. If the 
step direction changes sign, we prefer smaller moves. Secondly, 
if the function is exhibiting linear behavior on this move or 
thirdly, follows a sequence of best x values, a boost in the size 

of the step is desired. 

Linear behavior is measured using the formulas 
a = max(f0 - p2,0) 


a * (a + P*(A»,jC + - P))/(2a + P*ALFAC*(fj^ « P) ) 

ALFAC • 0.125 • (1/a + ALPAC) 

This fornula has the following effect. If f(? < p, then 

a » 1+1/ALFAC. If fg « P * fr, then a - 1. If fg » P*fi,» then 
u - 0.5 (One can think of increases in a as measuring more 
linear behavior). The factor for boosting the bounds is updated 
using the formula boost “ min ( (a*i‘l)*b^st, 10^0) . The for- 

mula for updating the bounds is 

(4) Bi » a*SI6N(6i’*'Bi,$i) boost*«i 

The sign of is used to store the direction of the last 
step. Mote that this formula includes the three effects men- 
tioned above. 

3. Selecting the Dimensions After Finding a Best Value of 
itf (x)fi After a Value that is not a Best 


The size of the move is based un differences between best x 

values using (4) with a ^ 1/4 and boost - 1. This choice was 
based on performance. A conservative measure is used because a 
new best was not immediately found. 

It is important to use the difference between best x values 

rather than the size of the individual moves. If adjacent best 
X values are smaller than moves taken, then there is nc benefit 

to taking large moves. On the other hand, if adjacent best x 

values are significantly farther apart than the size of moves 
taken, this is an indication that the large move can be allowed. 

4. Selecting the Dimensions After Givi ng Up and Retreatin g to a 
Best Value for iif(x)ii 


The logic for this case is the simplest of all. The dimensions 
of the box, after retreating, are reduced using (4) with 
a » 1/8 and boost 1/4. Some reduction is required given 
that a retreat took place. The values of 1/8 and 1/4 were 
chosen based on performance with the test set. The move Si 
is the last move that resulted in a new best value for 
iif(x)ff, or the value used after adjusting from a previous 

A# 

retreat. 
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5. Selection the Dinension After Obtaining a Value for itf (x) H Larger 


than th« Best Obt«li»«d Ptwleallv 

In this case it is in^rtant that the sise of the box should tend 
to decrease. If the sise decreases too nuch^ however, out first 
rule for retreating would take effect and waste the efCort ^e 
formulas for the else of the bound in this case are olven by 

a = {(l/4)*f^ + (l/2)*fj^)/(f^ 4 t^) 


(5) B. » a*B., if • 0 


« a* [SIGN ((1/4)*.^. + B^,5^) + (3/4)*«^J, ^ 0 


This formula has the following properties in each separate coor- 
dinate of the box. If the move gave a good decrease in f^ and 
the direction of the move is unchanged, then the dimencfon of the 
box stays the same. If £q is a good decrease in t^»e opposite 
direction, the box shrinks by 3/4. If fc significantly in- 
creases and the move is in the sane direction, the box shrinks 
by 1/2. Finally, if f^ significantly increases and the move 
is in the opposite direction, the box is reduced by 3/8. 

We emphasize here that it is not necessarily a bad move if the 
value of f.Q increases but the linear residual P is still small. 

4 . Some Test Results for the New Algorithm 

We developed pur algorithm using the test set [5] for uncon- 
strained nonlinear least squares problems. Raving this set of 
testing software was most appreciated by the authors. In this 
section, we give the results of our algorithm compared the 
MINPACK [1] codes for nonlinear least squares problems. 

There are characteristics of some test problems that make com- 
parisons hard to draw. For example certain problems, (8, 9, 10, 
15 (10 by 10), 16 (N by N) , N*10, 3C, 40) have more than one 
local minimum. The authors of [5] claimed that a package they 
tested, NLSQ2, failed on problem 10, second starting point. This 
problem has N « 3 variables and M » 16 data values. In fact, 
the result corresponds to a local minimum near X2 * X3 <■ ». The 
residual reported for this problem is the residual variance 
obtained using the average of the data values. For the class of 
algorithms being considered, we would not characterize any con- 
vergence to a local minimum as a failure. 
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For these reasons, we present results in Tables 1 and 3 for the 
B-K~l and NIMFACK codes for those problems that are directly 
comparable, in particular, the points that the two codes con- 
verge to are the same, and the convergence criteria are 
(effectively) the same. The output flags INFO/H-K-1 and INFO/ 
MINPACK are defined below. Their values indicate the reasons why 
each algorithm converged on each problem. 

Tables 2 and 4 show the results of the H-K-1 and MINPACK codes 
for those problems that have the same convergence criteria but 
with different solution points obtained. Note that the distance 
of the solution from the initial point is listed for both 
solutions. 

Reliability of software for nonlinear least squares has been 
discussed in the evaluation paper [6] . A reader may get the 
impression from the paper that efficiency is not an issue in 
evaluating code packages, only reliability. Obviously, effi- 
ciency of the algorithm is important. One of us (Krogh) has been 
concerned with applications involving multipoint boundary value 
problems with modest sized systems of ordinary differential 
equations. Each evaluation of the function and Jacobian matrix 
requires the integration of the system over a long trajectory, 
and the computation is expensive. In fact, this part of the 
process dominates the linear algebra involved in the Levenberg- 
Harquardt algorithm used to implement the Gauss-Newton method 
[7]. For problems where the functions are cheap to evaluate, the 
issue of efficiency in the linear algebra can become important. 
Of course, reliability is important. The software should not say 
that it has found a solution when it has not; it should not fail 
because an Intermediate point has a singular Jacobian matrix; it 
should satisfy all of its claimed tests for convergence, etc. 
It is necessary that a careful study of the output from a package 
be made so that fair conclusions are obtained when one makes 
comparisons. For example, all of the packages evaluated in [6] 
converge to local minima; none of them guarantee that a global 
minimum has been found. Thus, one cannot fault any algorithm 
tested for finding one local mininum versus another. (One might 
prefer the algorithm that moves to the nearest min.) Our view 
is that performance should be compared only for those problems 
that result in convergence to the same local min. 

We used constraints only for problem 10. (Results are also given 
without constraints.) Constraints used were: (a) nonnegativity 
on the variables and (b) restricting the argument of the expo- 
nential functions so that the value is within machine limits. 
These constraints seem natural for this problem. From the point 
of view of a user of nonlinear least squares software, one would 
probably solve the problem without constraints at the outset. 
Then after seeing the unsatisfactory results, constraints would 
be applied. After this the user would have more confidence that 
the desired minimum value had been obtained. 
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We now give a list of meanings for the output flag IMPO « The 
machine precision was n * 7.1 x on a CDC 6600. 

INFO/H“K-l Meaning 

2 Obtained a small value for iif(x)ii. (Tolerance is 

nl/2). 


3 Obtained a minimum for tif(x)ii. 

A# Ai» 

Tests satisfied are: 




(a) 

max(0,f£ - 

p2)l/2 < 

10-2*p 


and 

(b) 

life - fell 

< 

10“2*fg 


an? 

(c) 

life - fL'l 


fB*hl/2 

4 

No 

projection of 

residual 

into column space of 


Jacobian matrix, i.e., P > fc* (Usually means 
residual is pure noise.) 

5 liAxpil < n and retreat = .false. 

6 iiAXpii £ (10n)*iiXpil and retreat = .false. 


INFO/MIWPACK Meaning 

1 Obtained a minimum for iif(::)ii 


Test is: 

fc 1 U + 

2 llAXpll < (10n)*IIXBll 

• ^ 

3 INFO/MINPACK = 1 and 2 are both satisfied 

4 No projection of residual into column space of 
Jacobian matrix. 

5 Obtained a small value for iif(x)ii. (Tolerance 

is • 


The results listed in Tables 3 and 4 are only for the reader's 
general information. No conclusions should be drawn from the 
numbers. We plan to present algorithms that will solve the 
problems 6, 7, 13, 14 and 15 in Part II of this report (not yet 
written) . We expect the performance of our algorithm to improve 
most markedly for the nonsmall residual problems. 


Interpretation of the results of Tables 1-2, vith regard to 
efficiency, can vary. For exanple, problems where the Jacobian 
matrices are analytically available frequently have the property 
that evaluating f(x) and 3f/3x involve little more work than 

A# ^ MM 

evaluating f (x) alone. In many applications these quantities 

can be efficiently computed at the same time, and users often 
will do so. In this case, ' the amount of work for each problem 
would be proportional to the numbers NFEV in Tables 1-4. Under 
this measure for efficiency, our new algorithm is (almost) uni- 
formly doing a superior job since the function and Jacobian are 
required at the same point. At the other extreme, in applica- 
tions where the Jacobians are approximated by one-sided divided 
differences, the amount of work for 3f/3x will be about N 

M M 

times the work for evaluating f(x). For this comparison we com- 

M M 

pute WORK « NFEV -I- N*NJEV, summarized in Table 1 and 3. Using 
this measure for work penalizes our performance because of the 
fact that at each--— iteration we evaluate both f(x) and df/dx. 

MM MM 

Even under thes^ conditions we are generally performing in a 
Sv'perior manner with regard to efficiency. 

Th> problem listed under NPROB in the tables has either a blank, 
a ' or a ” by its number. These denote the standard starting 
point, [5] , 10 times it, or 100 times it, respectively. 


RJH : sdc : 7289A: 08/07/84X 
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ConqEMirable Results for small Residual problems 
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